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BAND GAPS IN 2D PHOTONIC CRYSTALS WITH 
HEXAGONAL SYMMETRY 


Bogdan LAZAR’, Paul STERIAN” 


Rezumat. Articolul expune, in detaliu, 0 procedurad matematica ce poate fi utilizata la 
calcularea  diagramelor de _ dispersie  corespunzdtoare  Cristalelor Fotonice 
bidimensionale cu simetrie hexagonala: In final,.sunt date citeva rezultate numerice 
pentru a dovedi validitatea metodei prezentate. 


Abstract.The article presents, in detail, a mathematical method useful for calculating 
dispersion diagrams corresponding to Photonic Crystals with hexagonal symmetry. In the 
end, a few numerical results are given to confirm the validity of the method. 
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1. Introduction 


Mathematical calculations and practical experiments show that a composite, 
formed by a repetitive succession of media with different dielectric permittivities, 
named also Photonic Crystal, possesses frequency gaps and as a result the 
electromagnetic fields, with speeds of oscillation inside those gaps, can not 
propagate through it [1], [2], [3], [4]. Therefore, photonic crystals can be defined 
as periodical media that have the property of forbidden frequency ranges, a 
radiation with the wavelength in their frequency gaps being unable to propagate 
inside these repetitive composites. 

The most usual and interesting type of photonic crystal, to date, is a dielectric 
material characterized by a cyclic electric permittivity that repents in space with a 
period comparable, as linear dimensions, with the wavelength of the radiation 
interacting with the dielectric. 

No simple formula, able to predict the size and positions of photonic crystals band 
gaps, exists [5], [6]. Unfortunately, when it comes to establishing the dispersion 
diagrams of this type of alternating structures, various articles present the results 
specifying that they have been obtained using a certain numerical method (for 
instance PWM — Plane Wave Method) implemented with a software conceived by 
the author, which if available is not well documented and written in a language 
you are not familiar with. For this reason, programs that calculate the structures of 
forbidden bands are hard to integrate in your own software, designed to study 
various properties of photonic crystals, and in conclusion, many people have to 
write their own piece of code able to calculate the dispersion diagrams, in other 
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words, to solve Maxwell Equations for a periodic dielectric medium in the 
frequency space. 

The purpose of the present paper is to start from these equations and finally get a 
mathematical set of expressions that can be easily implemented in software, 
especially Matlab, with the goal of obtaining dispersion diagrams for any 2D 
dielectric photonic crystal having hexagonal symmetry. The cases of 1D crystals 
[7] and 2D structures with square symmetry can easily be particularized from the 
hexagonal one. 


2. The wave equation in the frequency domain 


In classical physics, the propagation of electromagnetic waves in substance is 
studied using Maxwell Equations. The photonic crystals, being a repetitive 
succession of media, each of them extending in a volume many orders of 
magnitude greater than the dimensions of atoms, are perfectly suitable to be 
treated with these equations whose general form is: 


vxE- 2B (1) where: E=E(r,f) is the intensity of the 
a” electric field, B=B(r,t) the magnetic 
V-B=0, (2) induction, H=H(r,r) the intensity of 
A magnetic field, D=D(r,t) the electric 
VKH =~ + i(r.), (3) induction, j(r,t) the current density 
Ot ‘ : 
and p(r,t) the electric charge density. 
V-D= (r,t), at ) 


In Cartesian coordinates, the position vector r has the expression r=xe,+ye,+ze;, 
where exz,y are versors corresponding to the x, y, z spatial directions. 
The quantities D and H are, in general, for an arbitrary medium, complicated 
function of the following four variables: t, r, E and B: 

D=D(.r;E,B), H=H(+r;E,B). (5) 
However, for an entire group of substances, relations ( 5 ) turn into simple linear 
dependencies if the intensities of E and B are relatively small. Thus, 

D== €)¢,E, t--™= z : 

Ho Molly 
where: € is the electric permittivity of the medium, w- magnetic permeability, «, 
Ho - electric permittivity and magnetic permeability of vacuum respectively and ¢,, 
Hy - electric permittivity and magnetic permeability of the medium in respect to 
vacuum. The equations ( 1 ) - (4 ) can have an even simple form if none of the 
substances under consideration is magnetic, 
HM, =1, (7) 

and no density of electric charge or current exists, 


j(r,1)=0, p(r,r)=0. (8) 


(6) 
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Conditions (7 ) and (8 ) are met for the majority of dielectrics, at small intensities 
of electric and magnetic fields. Unfortunately, all simplifications end here because 
photonic crystals have position dependent electric permittivity in the form of a 
repetitive function of r: 
é,(r)=<¢,(r+A), where A is the period. (9) 
Therefore, substituting relations (6 ) - (9 ) into Maxwell Equations ( 1 ) - 
(4 ) and solving the system, two equivalent propagation equations: ( 10 ) and ( 11 
) and two conditions: {( 12 ),(-14 )tor {(.13.), (15 )}.can be written: 


1 1.0°?H 1 1 OE 
Vx VxH /=- i 1 Vx(VxE)|=-——_, 11 
(sa ce ar2 ( 0 ) &, 6} ( ) 2 ar ( ) 
V-H=0, (12) V-H=0, (13) 
V-[e,()E]=0, (14) V-[e,@)E]=0, (15) 


In the case of photonic crystals, of great interest is solving the two equations in 
the frequency domain. For that purpose H and E are considered being harmonic: 


(a) H(r,t)=H(r)e/™,  (b) E(r,t)=E(r Je! , (16) 
where @is the pulsation of the field: 
o=2af. (17) 


By replacing ( 16 )(a) in ( 10 ) or ( 16 )(b) in ( 11 ) the following atemporal 
expressions are obtained: 
2) 
ve{ ovat) |= 2H), (18) : (vx [vxE(r)) =? Er). (19) 
oF (r) c] é; (r) oe 

The equations ( 18 ) and ( 19 ) are useful for calculating the dispersion diagrams 
of various photonic crystals and implicitly for establishing their structures of 
frequency gaps. The equality ( 19 ), for instance, is general and can be solved for 
the full 3D case or simplifications in one or two dimensions. As stated in the 
beginning, this article deals only with the circumstance where the medium, taken 
into consideration, is two-dimensional, a case that splits in two branches. 


The first is the transversal electric possibility when E= E--e, that once replaced in 
( 19 ) transforms it in ( 20 ) (the e, versor and z index will be considered implicit). 


1 a” E(x, y) 0° E(x, y) oO 
weal are ay? a2 (x y) or modes (20) 


The relation ( 20 ) belongs to a category of equations that can be solved using the 
Bloch-Floquet theorem which (for the current situation) states that: if 1/e,(x,y) is a 
periodic function then: 


Ex, =e”) (x, y), (21) 
where g(x,y) is a repetitive function having the same period as 1/e,(x,y). 


In the particular situation of photonic crystals, ¢{x,y) is by definition periodic 
which also imply that 1/e,(x,y) is also cyclic with the same period as ¢,(x,y). 
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The second case is the transversal magnetic one, when H(r)=H-(x,y,z)-e,. This 
time, the equation ( 18 ) is used. First of all, the quantity V x|(1/<,(r))V x H(r)] needs 
to be evaluated. Thus, 


! __@|_ 1 oH,), @/ 1 GH,) 
vo arn))--2l 5S Peabo = CH 


Therefore, the following equation in H is obtained (where the z index and e, 
versor are considered implicit): 


af 1 \eH(xy) Of 1 .)eH(ey), 1) (Ay) Ay) ao” 
Zs) Ox Ase) oy wall ax? . oy a Hi(ey) ( = ) 
which can be used for calculating dispersion diagrams for TM modes. 


to 


3. The Fourier Transform for periodic electric permittivity 


The equation ( 20 ), with ¢, repetitive, can be solved by applying the Fourier 
Transform to both sides of the equality. In the case of photonic crystals with 
hexagonal symmetry the basic brick of the structure is rhombic, like in Fig. 1, and 
in consequence the periodicity of ¢,(x,y) can be mathematically written as in 
expression ( 24 ). 

The present paragraph calculates the Fourier Transform of ¢,,, 1/e, or in general, of 
any repetitive 2D function having a periodicity as that in Fig. 1. Once this 
transform is found, it can be used for solving equation ( 20 ). 


/e/ 
LOYOL 


Fig. 1. Photonic crystal with hexagonal symmetry. 


é(x, rae actal 140s Zhe tasing | (24) 


By definition, a square integrable function with n variables can be written as an 
integral sum in the following form [8]: 


1 t t t AUaens x. x, 
Pata) =F [f+ [a Gak ky, eee die dkig dls (25) 


kyy=-Oky»=-00 kk, = 080 
where p is the Fourier Transform of jf. Consequently, using ( 25 ), where f is 
replaced by e, both members of the equality ( 24 ) can be expanded as follows: 
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1 i ee2@ Voe.( yoo 
WG : Joa[os 2 ot, Vuy 


i(k,x+k,y) 1 +f 
dk dk, = ) 

oe =—00 k,,=—00 t =o ky 
( 26 ) 


which is satisfied if: 


3 k.=n an 
a | a 
i (e/344,) a3 a 
e 2 =1> (c, V3.nke l= 2nz > . re Qn’ n,, N» integers. (27) 
y Ba 
As a result, ¢(x,y) is constrained to have the decomposition: 
PR, ed 2 2 2 jn te | a 4 
a 2N,—Nn, 20 allt H 40 
E(x, y= n, —,—*=—- — le (28 ) 
oD ae DD Nee g | a al 
= 22 ms Ba »| 
and integrating over one 


By multiplying both members with © 
period (see Fig. 1), the expression ( 28 ) turns in: 


3 at 
2 3 j28{ mx 2 7 
| fel x, tee v3 dxdy = 
y 
i & 
V3 4 > > (29) 
2400 2 22 2ny-m 2a a 8 PE om peace = 

Sg —_——_ ee dxdy 

Mae Bee ad, | 
> 2m 2n.-n a “a 

Aa a’ iB a if 

where J is 
e27a 1274 
ae ce i\ V3 
a v3 
pS u 5 pt.my =m), Ny = My (30) 
Om > 0 for the remaining cases 
In conclusion 
Bx : 
7 3 - 58 | ma »| 
ax ye ‘3 dxdy ~ 


: 22 2m,—m, = 
1 5) 
a Wi a 
(31) 
N-1N-1 avs) Falmer "| 
a |< e * 


q=0 r=0 


Using (31 ), equations ( 20 ) and ( 23 ) can be solved (see the next paragraph) 
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4. The matrices that give the dispersion diagrams for TE, TM modes 


TE modes: Using (28 ), ¢,'(x,y) from ( 20 ) can be put in the form: 


24| (2n3—-n, ) 
jim— i y 
é, (x,y)= =D > Paint 

Ny =-0 Ny = 


(32) 


Also, E(x, y)= oilkertyy) (4, y) is expanded as follows: 


Tere 


PED! Re Pe : (33) 


, are the coefficients of the two series. 


E(x, y)= 


where p,, », and h 


my mM 


Thus, 0°E(x, y)/ax? + 6?E(x, y)/@y? can be written as: 


2 , F oo 2 me 2 ( an oh aes igri 
Blxy) | OPEC y) __ 2 y ay Iona rn, ( +E m | (1,42 Seren) i 1 : Al V3 } (34) 


ex? dy? ENB er v3 


and so, ( 20 ) transforms into: 


E & z 2a (|, 20 (2my=m) ake. [oma macmieotary| 
¥y p2. », yy lO cr Ph k, 4 mM ky } TB ie 


3 ‘sz = a 
Ny =—90 Ny =—20 Ny =—910 mM, =—910 


8 


(35) 
oe a 2a 2m ei 2m,-—m, ; 
ne) 
ao 3 
Res My 


Mee both members of (35 ): with (1 a? sin(z / 3) ede! a)nix+(2ms—mi)»/V3] and 
integrating with respect to y over the interval [ -ay3/4, av3/4] and x over [ y/¥3 - 4/2, 
y/V3 +a/2], the following equality is found: 

a we » ‘ _ 2a (2m, 1) sin[z(m, +n, —m; )] sin[z(m, +N, —m, )] _ 
Dy > y¥ “ey M,N Pn, sNy ( ‘ | u i, i 3 | 7 7 


ny =—00 ny =—00 my =—00 a a(m, +n, — mt) am my +N —m}) 


om, +N my ing +N, 


Ce, 00 am = m}) am, -m’) 


Onn, my ny my 


oO = x .. sin[z(m, - m')|sin[z(m, - m',)} 


(36) 
For an arbitrary (mj,m) pair, the majority of terms, situated left and right in respect 
to the equality sign, will disappear and (36 ) simplifies to: 


00 00 2 2 2 
20 20 (2m, —m,) @ 
. ) Fn JN Prnt—m,,my—-m, G + 4 m o{s + B = 2 Pont > (37) 


a 
Mmy=—% Mm, =—00 


for TE modes 
which represents a system of equations that, if solved, gives a set of eigen 
frequencies, wx. 
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In practice, the m indexes will be taken: m,,m{,m,,m, © [-M, M] where M is a 
positive integer. For each of the (2M+/) x (2M+1) values of (m,m5) an equation 
like ( 37 ) exists where the coefficients p,._,, m,-m, have indexes that vary in the 
interval [-2M, 2M] and can be calculated (see (31 )) with the formula: 


. 2a (2n, n) 

1 4M 4M : av3 =i | mas a| 

Pry =o Df | eh Th a) ne (38) 
(4M +1) soy | | 4rd "2(4M +1) °° 


where the following notations have been made: n, =m|—m, n, =m, —m,. 


In conclusion (37 ) with m,,m;,m,,m, © [-M, M] is a system in the form: 


Cc 


2) 2 
@ 
S-h=“—h or [s- 0, (39) 


where S is a square matrix having (2M+1)* elements of the type: 


2 2 
20 2m \2m,—m 
Pin! —m, ,m’,—m, (: + aa) fe, a 2 Ens | | ( 40) 


that can be calculated for any given m,,mj,m,m},k,,k, . 


of unknown 


values. It can be noticed that ( 39 ) is satisfied, independently of h, if 
det(S <pe)= Owhich gives (2M+1)* possible ax. 


h is a column matrix possessing (2M+/) x (2M+1/) elements, h 


™M,INy ? 


Therefore, for any given (k,, ky), (2M+/ y° values for @ are found and, in this way, 
the dispersion diagram w=a(k,, ky) is obtained. As can be noticed, a single pair of 
given (k,, ky) require solving a system of (2M+/ Y equations where for good 
precisions M have to be increased till no difference is observed between @= a(k,, 
k,) calculated with (2M+/ yY and with (20(M+/)+1 ‘ equations. 


TM modes: The expression ( 37 ) is valid only for the TE modes. For 
finding its equivalent corresponding to the TM situation, the equation ( 23 ) have 
to be utilized as starting point. Using the same deductions as in the case of TE 
modes, the following expressions can be successively written: 


1 0°H(x,y) 0°H(x,y) _ 
ax? dy? 


OT eee ee ee ee en 2 On (2 ei ) 2 jee en | 
= 2/3 y os > » Hains Deuce [, m ic + ic l | & 
a n = 


| =—-00 Np =—00 MM, =—20 My =—00 


(41) 
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°/ 1 Lo a F : jae. 


éx dy e,(xy)) ay 


2 co co on 128m Maat ) 2 on iere22{ mae 2 My) 
j— ‘ hi jk, +— ‘ 


Mm =-2 My = 


Edm i : ii apte), § 7 : ‘ [2mm boston] 
Se a a 3 a V3 

(42) 
The decompositions obtained from (41 ) and ( 42 ) together with ( 32 ) and ( 33 ) 
(where E(x,y) is replaced by A(x,y)) are introduced in ( 23 ). After simplifying by 
eK", multiplying by (/ a’ sin(a/ gy) /07! alixems—miv/B] and integrating in respect 
to y over the interval [ -cv3/4, 23/4 | and to x over [ y/V3 -a/2, y/V3 +/2], the following 
equality results: 


$ > > > —— atan 2 m( bem) 2 Grate, «28 set om) of, tee | 


Ny =—00 My =—00 mM, =—90 My =—00 


_sinie(my +n, —m!')]sin[z(m, +n, —m))| e o y x h sin[z(m, =m{)]sin[z(m —m | 
alm +ny—mji) —x(m,+ny—m) ec? m=20 m,=—00 mins alm, —my) (im, - m4) 


Sn any ml 8 ny eng mh Sin, ant myn) 


(43 ) 
For an arbitrary pair (m},m), the majority of terms from the left and right of the 
equal sign disappear, ( 43 ) turning in: 
ay ; my k, 7m] Qn 2m), —m,)—(m =m, Ze (2m, = é 
@ o 


=\ & - a a B °C ee , 
D3, , ea ees x 2 alee Air die 
7 ( In ) ( 2a (2m, — ) e 
+| k, 4 m, leet 
a om, Wi 
for TM modes 
(44) 


where the same explanations as given for equation ( 37 ), corresponding to the TE 
situation, remain valid. 

As already explained, by solving ( 37 ) and ( 44 ) the dispersion diagrams, 
®=(k,, ky), are obtained. It can be shown that, for a photonic crystal with 
hexagonal symmetry (whose basic cell is described by the vectors a; and a; as in 
Fig. 2), a reciprocal cell, defined by by, bz (see Fig. 2), exists in the spatial 
frequencies domain (dual space). Therefore, the periodicity e(r)=e(r+A) has a pair 
in the dual domain, o(k)=@(k+B), and in consequence it is enough to compute @ 
for k inside just one elementary cell B. More than that, if e(r) has some symmetry 
inside the photonic crystal cell, then also m(k) has symmetries inside B and this 
property further reduces the range of k for which @ have to be evaluated. In the 
particular cases of the crystals given as examples, in the following paragraph, it is 
sufficient to calculate the dispersion diagrams just for values of k lying inside the 
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triangular domain LXM (see Fig. 2), called irreducible Brillouin zone. More, 
numerical calculations show that the worst scenario, with the smallest band gaps, 
happens for wavenumbers k along the contour LXML, and for this reason, 
diagrams @=(k) will not be represented for the entire surface of the triangle 
LXM but just for the contour LXML. 


Lae Qn 1 Qn 2 
a, =ae,; a, =al)—e,+——ey b, =—|e,-—~e,| b, =—-—=e, 
1 x 2 c ) 2 y 1 al Wi 2 a i y 
A=l,a, +/,a, where |,,/, integers B=n,b, +nyb, where n,n, integers 
L-> (0,0) 
b, ‘ 
mx 1 
‘xX X —| 0, —= 
/ x =) 
\ ij 
: l 2n 2x 
\ M >| —,—= 


Fig. 2. The vectors that describe: normal space (left) and the dual space (right). LXML is the path 
alongh which the dispersion diagrams will be graphed. 
5. Dispersion diagrams and band gaps. Numerical results. 

Using ( 37 ) and ( 44 ) the dispersion diagrams for the TE and TM modes, 
corresponding to a few particular geometries of photonic crystals with hexagonal 
symmetry, will be calculated. A total of three configurations are studied. For each 
situation, the entry parameters are given in the description of the case, beneath the 
figure or diagram. The significations of these parameters are: 


E,q inside thespecific geometric element (circle, triangle, hexagon) 


5) 


G)eyonyo=4 
E,» in therest of theelementary cell (called backroung) 
(b) f = the filling factor defined as the fraction between the surface of the 
geometric element and that of the entire cell. (c) fina= maximum filling factor 
attainable. Also, for each case, some specific parameters as r (the radius of the 
circular element in Fig. 3) and d (the length of the triangle or hexagon edges, Fig. 
6, Fig. 9) are given. Another important parameter is N x N = (2M+1) x (2M+1) 
(see the explanations for the equation ( 37 )) that represents the number of 
discretisation elements in which the basic cell of the crystal is divided. 


Regarding the diagrams in Fig. 4, Fig. 5, Fig. 7, Fig. 8, Fig. 10, Fig. 11 a few 
explanations have to be given: (1) A scaled frequency, wa/2mc, was represented 
on the y axis in order the dispersion diagrams could be read at any frequency and 
any elementary cell size. Thus, supposing that wa/2mc=0.4 (see one of the 
diagrams) which is equivalent to fa/c=0.4 (w=2nf) or a/=0.4, and knowing the 
crystal cell edge length, a=0.6 wm, then A=1.5 wm which corresponds to 
f-66 THz. (2) Both transversal electric and magnetic diagrams have been 
represented on the same figure for each crystal. 
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The curves with dotted line correspond to TM modes and the ones with solid line 
to TE modes. Also, for clarity, the TM and TE gaps have been marked using color 
bands and in the case of total forbidden frequency zones (for both modes) the gap 
was hatched with oblique lines. As can be seen from diagrams, various band 
structures are obtained when €,a, &,» and the geometrical figure shape, inside the 
crystal cell, are varied. 


34 
A 
2nr* a y | 
p= ; 25 k / 
a'v3 20 / y 
/ 
a a y 
Tmax =D > 15 é | 
2 a y 
10 y 
f, = a . d y 
is p <_ | 
1 5 10 15 20 25 30 35 40 45 50 
Fig. 3. Periodic element defined by: r=0.2a; NXN=3333 (two cases of (€,,, €;,) are considered). 


: wa/2nc : wa/2nc 


0.8" ee | 


L X M L 


Fig. 4. r=0.2a3 €2=5; &y=1; NxN=33 33. Fig. 5. r=0.2a; €,=12; ¢4=1; NxN=3333. 
38 , 
, 
30} a 
= F 
a 20 y 
Lr, = a 3 y 
max circumscrbed V3 2 10| d 
ding = 4 5 Fives =! 1 y. y : 
1 10 20 30 40 50 56 


Fig. 6. Periodic element defined by: r..=0.42a; NxN=37%*37 (two cases of (a, &:») are considered). 
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@a/21c 
0.8; - 


SZ 


Fig. 7. Fec=0.2243 &=12; &4=1; NXN=37%37. Fig. 8. rec=0.2243 €=1; &5=12; NxN=37%37. 


38 
6d? 30} 
f =i 

20} 

Tmax circumscrbed = S > 
10} 

a 2 

WV == om SS 
NN 3 ee 1 
{ 


Fig. 9. Periodic element defined by: r..=0.27a; NxN=37%*37 (two cases of (€;a, &:») are considered). 
@a/2nc 


0.8 
0.7 


0.6 


0.8 


0.5 


0.4; 


0.3 


0.2} 


Tec=0.2743 Eyg=123 &p=1; NxN=37X37. 


M L 


0.7 


0.6 


wa/2mc 


O6 | 


0.45 
0.3 } 
0.2 | 


0.1} 


M L 


Fig. 11. re=0.274; &a=1; €p=12; NXN=37*37. 
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Conclusions 


Systems ( 37 ) and ( 44 ) can be used for obtaining dispersion diagrams and 
implicitly band gaps for a variety of dielectric two-dimensional photonic crystals 
with hexagonal symmetry. Both systems are in a form that can be easily 
implemented in software, especially in Matlab where, due to the richness of the 
already existing subroutines, just a few program loops need to be written for 
computing the coefficients.in ( 37.) and ( 44 ) with which a square matrix is 
generated and finally the eigenvalues of it are extracted using a general function 
already available in Matlab. Each set of wm eigenvalues corresponds to a 
wavenumber, k=k,ex+k,ey, that can be chosen to vary along an arbitrary path or in 
a given domain. In practice, due to symmetry reasons, it is enough to take k along 
LXML path (see Fig. 2). 

As regarding the numerical examples, they are given just for demonstrative 
purposes, in order to show the correctness of the formula written in the current 
paper. For a thorough investigation, on how the size of certain band gaps are 
affected by the contrast between ¢€,, €,, and other parameters of the crystal cell, 
many diagrams have to be computed while a single entry value is varied in a 
certain range of interest. 
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